Classical simulation of infinite-size quantum lattice systems in two spatial dimensions 
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We present an algorithm to simulate two-dimensional quantum lattice systems in the thermodynamic limit. 
Our approach builds on the projected entangled-pair state algorithm for finite lattice systems [F. Verstraete 
and J.L Cirac, cond-mat/0407066] and the infinite time-evolving block decimation algorithm for infinite one- 
dimensional lattice systems [G. Vidal, Phys. Rev. Lett. 98, 070201 (2007)]. The present algorithm allows for 
the computation of the ground state and the simulation of time evolution in infinite two-dimensional systems that 
are invariant under translations. We demonstrate its performance by obtaining the ground state of the quantum 
Ising model and analysing its second order quantum phase transition. 
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Strongly interacting quantum many-body systems are of 
central importance in several areas of science and technol- 
ogy, including condensed matter and high-energy physics, 
quantum chemistry, quantum computation and nanotechnol- 
ogy. From a theoretical perspective, the study of such systems 
often poses a great computational challenge. Despite the ex- 
istence of well-stablished numerical techniques, such as exact 
diagonalization, quantum monte carlo [1], the density matrix 
renormalization group [2] or series expansion [3] to mention 
some, a large class of two-dimensional lattice models involv- 
ing frustrated spins or fermions remain unsolved. 

Fresh ideas from quantum information have recently led 
to a series of new simulation algorithms based on an effi- 
cient representation of the lattice many-body wave-function 
through a tensor network. This is a network made of small 
tensors interconnected according to a pattern that reproduces 
the structure of entanglement in the system. Thus, a matrix 
product state (MPS) [ ], a tensor network already implicit 
in the density matrix renormalization group, is used in the 
time-evolving block decimation (TEBD) algorithm to sim- 
ulate time evolution in one-dimensional lattice systems [5], 
whereas a tensor product state [6] or projected entangled-pair 
state (PEPS) [ ] is the basis to simulate two-dimensional lat- 
tice systems. In turn, the multi-scale entanglement renormal- 
ization ansatz accuratedly describes critical and topologically 
ordered systems [8]. 

In this work we explain how to modify the PEPS algorithm 
of Ref . [7] to simulate two-dimensional lattice systems in the 
thermodynamic limit. By addressing an infinite system di- 
rectly, the infinite PEPS (iPEPS) algorithm can analyse bulk 
properties without dealing with boundary effects or finite- size 
corrections. This is achieved by generalizing, to two dimen- 
sions, the basic ideas underlying the infinite TEBD (iTEBD) 
[9]. Namely, we exploit translational invariance (i) to obtain 
a very compact PEPS description with only two independent 
tensors and (ii) to simulate time evolution by just updating 
these two tensors. We describe the essential new ingredients 
of the iPEPS algorithm, which is based on numerically solv- 
ing a transfer matrix problem with an MPS. We then use it 
to compute the ground state of the quantum Ising model with 
transverse magnetic field, evaluate local observables, identify 



the critical point of its second order quantum phase transition 
and estimate the critical exponent (3. 

We point out that the algorithms of Ref. [ ] have already 
addressed the computation of the ground state of infinite two- 
dimensional systems, by analysing an infinite transfer matrix 
problem with a MPS. A major difference in our approach is 
how this is handled. Instead of DMRG techniques (which 
consider an increasingly large chain with a finite MPS), we 
use the iTEBD algorithm [ ], based on a power method that 
uses an infinite MPS (iMPS) from the start. This seems to 
significantly improve the results reported in Ref. [6]. 

Finite PEPS algorithm. — We start by recalling some ba- 
sic facts of the PEPS algorithm for a finite system [7]. Con- 
sider a two-dimensional lattice C where each site, labeled by 
a vector r, is represented by a Hilbert space = C d of 
finite dimension d, so that the Hilbert space of the lattice is 
Vc = ®fec V^. For concreteness, we address the case of 
a square lattice, with N x N sites labelled by a pair of inte- 
gers r = (x,y), x, y = 1, • • • , N. [However, the key ingre- 
dients of the algorithm for infinite systems to be considered 
here are still valid for any type of regular lattice.] The model 
is further characterized by a Hamiltonian H = ^ h\ rir2 ^ 

that decomposes as a sum of terms h\ rir2 ^ involving pairs of 
nearest neighbor sites r\,ri G C. A pure state G Vc of 
the lattice is represented by a PEPS, namely a set of N x N 
tensors {A^}f? e c> interconnected into a network V that fol- 
lows the pattern of C (Figs, l.i and l.ii). Tensor A^ dlr 
is made of complex numbers labelled by one physical index 
s and four bond indices u, d, I and r. The physical index 
runs over a basis {|s)}s=i,— ,d of V^, whereas each bond 
index takes D values and connects the tensor with the ten- 
sors in nearest neighbor sites. Thus, is written in terms 
of 0(N 2 D A d) parameters, from which the d N complex am- 
plitudes (5(1,1)5(1,2) • • ' 5(at,at)|^) can be recovered by fixing 
the physical index of each tensor in V and by contracting 
all the bond indices. 

Given a PEPS for some state |^o) G Vc, the algorithm of 
Ref. [ ] allows to perform e.g. the following two tasks: (i) 
computation of expected values (\I>o|0|^o) f° r a ^ oca ^ °P era_ 
tor O, such as the energy, a local order parameter or two-point 
correlators, and (ii) update of the PEPS after a gate g^ lT2 ^ has 
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Figure 1 : (color online) Diagramatic representations of (i) a PEPS 
tensor A su dw with one physical index s and four inner indices u,d, I 
and r; (ii) local detail of the tensor network V for an iPEPS. Copies 
of tensors A and B are connected through four types of links; (Hi) 
reduced tensor a of Eq. (2); and (iv) local detail of the tensor network 
S. 



been applied on two nearest sites r\ , r*2 G £. The second task 
can be used to simulate an evolution according to Hamiltonian 
H, both in real time and in imaginary time, 
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in the sense of obtaining a new PEPS representation that ap- 
proximates the states |\E^) and |^ r ). This is achieved by 
breaking the evolution operators e~ lHt and e~ Hr into a se- 
quence of local gates, using a Suzuki-Trotter expansion [10], 
and by updating the PEPS after applying each of these gates. 
In particular, one can approximate the ground state of Hamil- 
tonian H through simulating an evolution in imaginary time 
for large time r, starting from a product state |^o) (for which 
a PEPS can be trivially constructed). 

Let £ denote the network made by the TV x TV reduced ten- 
sors (Figs, l.iii and l.vi), 
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where u represents the double bond index (u,u') and the 
physical index s has been contracted, and let f\ , r*2 G C de- 
note two nearest neighbor sites. Then the environment £^ ri ^ 
for these two sites is the network obtained by removing the 
tensors and o\ r ^ from £ . By "contracting a tensor net- 
work" we mean "summing over all the indices that connect 
any two tensors of the network" . It turns out that both (i) the 
computation of an expected value (^|Ot ri ' r2 ] |\fr) and (ii) the 
update of the PEPS after a gate g^ 1 ^ can be achieved by con- 
tracting £[ r i> r 2]. However, the cost of this contraction grows 
exponentially with TV. The core of the PEPS algorithm [7] is 
an approximate, efficient (quadratic in TV) scheme to contract 
£[ri,r 2 ] ? based on MPS simulation techniques. 

Infinite PEPS algorithm. — In order to consider the limit 
of an infinite lattice, TV —> oo, where both and H are in- 




L n r 2 ] 



(ii.a) jrthh) 
C D C D C D 




Figure 2: The environment c^ 1 for a link of type r is first approx- 
imated by an infinite strip J 7 ^ 1 ' 7 ^ and then by a six-tensor network 
g[ri,r 2 ]^ jhgsg reductions can be performed according to either a 
vertical/horizontal scheme (ii.a) or a diagonal scheme (ii.b). Ten- 
sors A, A*, B and B* are not part of the environment. 



variant under shifts by one lattice site, we need to understand 
how to efficienlty contract an infinite environment £ \- ri ^ 2 \ 
Translational invariance allows us to represent the state 
in terms of only two tensors A and B that are repeated indefi- 
nitely in V (Fig. 1), 

A l(x,x+2y)] = ^ A [{x,x+2y+l)] = ^ x y ^ ^ (3) 

so that the iPEPS depends on just 0(D 4 d) coefficients. No- 
tice that £ t ri ' r2 ] is also made of infinitely many copies of just 
two reduced tensors a and b, defined in terms of A and B ac- 
cording to Eq. (2). Then its contraction is achieved in two 
steps, as illustrated in Fig. (2): first we approximate £\ ri ^ r " 2 \ 
with an infinite strip jF[ ri ' r2 ]; then we approximate J 7 ^ 1 ^ 
with a small set of tensors Q^ ri ^ = {Gi, • • • , G$}. This can 
be achieved using a vertical/horizontal scheme or a diagonal 
scheme as illustrated in Figs. 2-3. 

The first step considers a transfer matrix R consisting of an 
infinite strip of reduced tensors a and b (Fig. 3). R can be 
regarded as a linear operator acting on an infinite chain where 
each site is described by a vector space of dimension D 2 . Let 
|$) denote the dominant eigenvector of R — that is, the eigen- 
vector of R, R\&) = A|$), with the eigenvalue A of largest 
absolute value. Here we assume that the dominant eigenvec- 
tor is unique [11]. By construction R is invariant under shifts 
by two sites of the infinite chain - and so is |<I>). We use 
an iMPS, characterized by just two tensors {C, D} and with 
Schmidt rank to represent an approximation of |$). We 
obtain this iMPS by simulating (repeated) multiplication by 
R on an initial state |$q) w i m me iTEBD algorithm [ ] and 
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Figure 4: Transverse magnetization m x and energy per site e as a 
function of the transverse magnetic field h. The continuous line 
shows series expansion results (to 26th and 16th order in perturba- 
tion theory) for h smaller and larger than h c ~ 3.044 [14]. Increas- 
ing D leads to a lower energy per site e. For instance, at h = 3.1, 
e(D = 2) « -1.6417 and e(D = 3) « -1.6423. 



Figure 3: Transfer matrices R (i.a) and S (i.b) for the verti- 
cal/horizontal contraction scheme. Multiplication of an iMPS by R 
using the iTEBD algorithm comes at a computational time that scales 
as 0(x 3 D 6 -\~x 2 D 8 d) (similar costs apply to multiplying by S). This 
cost is lower in diagonal contraction scheme (ii.a) and (ii.b), namely 
0(x 3 D 4 + x 2 D Q d), but a slightly larger x is required in order to 
retain the same accuracy. 



by using the fact that 
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The iMPS for |<E>) accounts for an infinite half plane of the 
environment S^ ri,r2 \ Similarly, we use another iMPS with 
tensors {C, D'} to encode the left dominant eigenvector (<I>'| 
of R, ($>'\R = A( < l )/ |, which also accounts for an infinite half 
plane. Then jF[ ri ' r2 ] is built from these two iMPS and a strip 
of reduced tensors a and b. 

In the second step, a transfer matrix S is defined in terms of 
the tensors {a, 6, C, D, C , D'} (Fig. 3). S can be regarded as 
a linear operator acting on three sites with local vector space 
dimensions x> D 2 and x- Again, its dominant eigenvector 
\Q), encoded in a three-legged tensor X, is computed from an 
initial state |Qq) using the fact that 
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Let X' be the tensor for the left dominant eigenvector (Q'\ 
of S. Then Q^- rir2 ^ is a (circular) MPS consisting of the six 
tensors {C,D,C',D',X,X'}. 

Simulation of time evolution. — We decompose the 
Hamiltonian as H = H r + Hd + Hi + H U9 where the operator 
Hr = E(ri r 2 ) r h^- rir2 ^ collects all interactions along r-links 
(and similarly for d-,l- and w-links), and consider a Suzuki- 
Trotter expansion of the time-evolution operator e~ %Ht of 
Eq. (1) in terms of operators e - ^ r( ^, e ~^ d5 \ e~ lHlSt and 



e -zH u 5t^ w h ere §f [ s some small time step. Each of these op- 
erators further decomposes into a product of identical two-site 
unitary gates g = e~ zh6t acting on all pairs of sites connected 
by a link of the proper type. For instance, for links of type r 
we have 
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Without loss of generality, we need to address only the up- 
date of tensors A and B after applying e - lH ^ 5t to Let 
us assume that the gate g is applied on just one of the r-links. 
In that case, in order to update the iPEPS we would (i) com- 
pute the environment for that specific r-link following Figs. 
2-3, and (ii) determine the optimal new tensors A! and B' for 
the link, using the optimization techniques of [7] (sweeping 
over just the two sites involved). We notice, however, that the 
above A' and B' already define an iPEPS for e - iHr5 \^) - 
that is, with gates g acting on all r-links. Indeed, this follows 
from translation invariance and the fact that a unitary gate g 
on a given r-link does not affect the environment on any other 
r-link. In other words, the update of tensors A and B on each 
r-link is identical and need only be performed once. 

The above argumentation fails for an evolution e~ Hr in 
imaginary time, since the gate g' = e~ h6r is no longer uni- 
tary. In this case, a gate applied on, say, an r-link modifies the 
environment on the rest of r-links. Nevertheless, as in one- 
dimensional systems [9], the same algorithm can still be used 
to find the ground state of the system through imaginary-time 
evolution, provided that a sufficiently small St (leading to al- 
most unperturbed environments) is used at the last stages of 
the simulation. 

Quantum phase transition. — To demonstrate the perfor- 
mance of the iPEPS algorithm, we have simulated an evolu- 
tion in imaginary time to obtain the ground state |^a) of the 
quantum Ising model with transverse magnetic field, 
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Figure 5: (color online) Magnetization m z (X) as a function of the 
transverse magnetic field A. Dashed lines are a guide tp the eye. 
We have used the diagonal scheme for (D,x) — (2,20), (3,25) 
and (4, 35) [12] (the vertical/horizontal scheme leads to comparable 
results with slightly smaller The inset shows a log plot of m z 
versus |A — A c |, including our estimate of A c and f3. The continuous 
line shows the linear fit. 
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Figure 6: (color online) Two-point correlator S zz (l) near the critical 
point, A = 3.05. For nearest neighbors, the correlator quickly con- 
verges as a function of D, whereas for long distances we expect to 
see convergence for larger values of D. 



Then we have computed the energy per site e and the trans- 
verse and parallel magnetizations m x and m z (Figs. 4-5), 

m x (X) = (^a|^|^a), m z (X) = (V x \a z \V x ), (8) 
and the two point correlator S zz (l) (Fig. 6) 

S zz (l) = <*A|4^f + ^ ] l*A) - (m z ) 2 . (9) 

Comparison with series expansion results of Ref. [14] 
shows remarkable agreement for all local observables on both 
sides of the critical point, which Monte Carlo calculations [13] 





QMC Ref. [ ] 


D=2 iPEPS 


D=3 iPEPS 


D=3 VDMA Ref. [ ] 


A c 


3.044 


3.10 


3.06 


3.2 


P 


0.327 


0.346 


0.332 





Table I: Critical point and exponent /3 as a function of D. 

indicate to be at magnetic field \mc ~ 3.044. We also obtain 
accurate estimates of the critical magnetic field A c and critical 
exponent (3, which for D = 2 and D = 3 agree with Monte 
Carlo results within 5.8% and 1.5% respectively. 

In conclusion, we have presented an algorithm to simulate 
infinite two-dimensional lattice systems. We have tested its 
performance in the context of the quantum Ising model, where 
our results can compete quantitatively with those obtained us- 
ing long-established methods, such as quantum Monte Carlo 
[1] or perturbative series expansions [ ]. TheiPEPS algorithm 
can now be applied to address problems beyond the reach of 
quantum Monte Carlo (since it has no sign problem) and se- 
ries expansion methods (since it does not rely on an expansion 
around an exactly solvable model). Thus we expect it to be- 
come a useful new tool in the study of strongly interacting 
lattice models. 
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